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SUMMARY 


The method of lines (MOL) has been used to obtain solutions to axisymmetfic and 
two-dimensional inviscid flow over smooth blunt bodies. Comparisons with experimental 
data and the results of other computational methods have demonstrated that very accurate 
solutions can be obtained by using relatively few lines with this approach. The method is 
efficient; typical converged solutions have been obtained in 100 to 150 total integrations. 
The method of lines is semidiscrete and has relatively low core storage requirements as 
compared with fully discrete methods (such as time -dependent techniques) since very 
little data are stored across the shock layer. This latter feature is very attractive for 
three-dimensional problems because it enables computer storage requirements to be 
reduced by approximately an order of magnitude. 

The disadvantage of the method of lines is that roundoff errors increase exponen- 
tially with number of lines; thus, the total number of lines that can be used for a particular 
problem is limited. In the present study it was found that 9 lines was a practical upper 
limit for two-dimensional and axisymmetric problems. This condition limits application 
of the method to smooth body geometries where relatively few lines would be adequate to 
describe changes in the flow variables around the body. 

Extension of the method to three dimensions is conceptually straightforward; how- 
ever, three-dimensional applications would also be limited to smooth body geometries 
although not necessarily to a total of 9 lines. 

INTRODUCTION 

Over the past two decades, much effort has been devoted to solving the inviscid 
flow field over blvuit bodies. Two basic approaches have been used: (1) an inverse 
approach - where the shock wave is given and the body shape and details of the flow 
within the shock layer are unknown, and (2) a direct approach - where the body is given 
and all details of the flow within the shock layer are unknown. (See ref. 1.) 

For the inverse problem, the steady flow equations are usually solved by marching 
inward from the known shock wave to obtain the solution for the body. (See refs. 2 to 8.) 
Since the steady flow equations are of a mixed elliptic-hyperbolic type, this initial value 
problem is ill posed and inherent instabilities are present which complicate the solution 
procedure. (See ref. 9.) The inverse problem can be made direct by iteration; however. 



it is not always clear how a solution which gives an approximation to a desired body 
shape is to be improved (ref. 1) and thus solutions for general body shapes may be diffi- 
cult to obtain. The inverse approach does have two attractive features; (1) solutions 
are usually very fast; and (2) relatively low computer storage is required, since it is 
unnecessary to store all flow field data within the shock layer. 

More recently, time -dependent finite-difference techniques have been used suc- 
cessfully to solve the more difficult direct problem. (See refs, 10 to 13.) The tech- 
niques take advantage of the hyperbolic nature of the unsteady flow equations by march- 
ing in time from approximate initial conditions to the desired steady-state solution 
(asymptotically). This initial value problem is well posed and thus eliminates the insta- 
bilities that arise with the inverse approach. In addition, the time -dependent technique 
can be used to compute the flow field around complicated geometries and the extension to 
three dimensions is straightforward. (See refs. 14 and 15.) However, these techniques 
require relatively large computir^ times and storage, since even when shock fitting is 
used (refs, 11 and 12), the entire shock layer must be discretized. 

Telenin (refs. 16 and 17) has proposed an alternate method of solving the direct 
blunt body problem called the "method of straight lines" which combines many of the 
advantages of the inverse approach with the ability to treat more general geometric 
shapes. This same method has been used to solve nonlinear conical flows in refer- 
ences 18, 19, and 20 wherein it was called the method of lines (MOL). 

The purpose of the present paper is to explore the application of the method of 
lines for computation of the axisymmetric and two-dimensional flow over blunt-nosed 
bodies immersed in an oncoming supersonic stream at an angle of attack of 0°. 

SYMBOLS 

A,B,C parameters in equation (30) for body shape 
a speed of sound, a/v^ 

B^ body bluntness parameter (see eq. (26)) 

Cp specific heat at constant pressure 

Cy specific heat at constant volume 

g,Gj^,G 2 ,. . .,Gg parameters defined by equation (12) 
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total enthalpy, h/v 
static enthalpy, h/v 
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indicator, j = 0 for 2-D flow, j = 1 for axisymmetric flow 


Mach number 


line number 


NETA number of integration steps from shock to body 


number of lines 
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OO 




pressure, p/p^V^ 
radius of curvature of nose, m 
polar coordinates (see fig. 1) 
time, 

dimensional time, sec 

velocity component normal to body, 

free-stream velocity, = 1 

dimensional free-stream velocity, m/sec 

velocity component normal to shock wave, ^v^) 

velocity components in polar coordinate system (see fig. 1), 

Cartesian coordinates (see fig. 3), y/R^, 

location of pole of coordinate system (see fig. 3), 



shock wave angle 


r 

’'o 

Ar 

A?7 
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ratio of specific heats, 
constant in equation (19b) 
local shock-layer thickness (eq. (8)) 
integration step size in -rj -direction 
step size in ^-direction 

( fi 7 

typically 10" or 10" 

convergence criterion 

criterion for modified iteration procedure 


angle between tangent to body and line r = Constant 
^ angle between tangent to shock wave and line r = Constant 

S> 

^r’^0 derivative of rj with respect to r and <p (eqs. (7a) and (7b)) 

|,7j transformed coordinates defined by equation (6) 

derivatives of ^ with respect to r and (p (eqs. (7a) and (7b)) 
p density, p/p^ 

Superscripts: 

(“) dimensional quantity 

(~) predicted quantity 

( )' derivative with respect to r? in equations (19) 
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Subscripts: 


b body 

s shock 

oo free stream 

n nth line 

max maximum 

METHOD 


The problem is to compute the inviscid two-dimensional or axisymmetric flow over 
a blunt-nosed body immersed in an oncoming supersonic stream. Solutions are to be 
obtained throughout the subsonic region and downstream into the supersonic region 
(fig- !)• 

Flow Field Equations 

The partial differential equations governing the two-dimensional or axisymmetric 
flow of an inviscid, nonconducting fluid in polar coordinates (fig. 1) are 
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Dt r P dr 
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denotes the substantial derivative and j = 0 for two-dimensional flow and j = 1 for 
axisymmetric flow. 

Noting that entropy is constant along streamlines leads to the following relation: 


i_ DE = (4 

a2 Dt Dt 

which can be used to replace the density derivative in equation (1) to obtain the following 
relationship: 


Dp 

Dt 


+ 


\ar r 30 r 



= 0 


(5) 


Equations (2), (3), and (5) contain derivatives of p, v^, and v^ only. 

Computational Domain 

To facilitate integration of this system of partial differential equations, the region 
(fig. 1) bounded by the stagnation streamline, the bow shock wave, a line ^ 

TllciX 

(located in the supersonic region), and the body surface is mapped into a rectangular 
domain (fig. 2) using the transformation equations 


^ = 


^max 




r - r^(0) 

rg(0) - rj^(0) ^ 


( 6 ) 


From these transformation equations, a set of transformation operators can be defined as 


i_ = 7] i- + f l- = JLl. 

3r ^ ^ 3| Ar 37/ 


I- = r] — + f — 

30 d7] 


(7a) 

(7b) 


where 


Ar = r^(0) - Tjj(0) 


(8) 
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which will transform the governing differential equations (2), (3), and (5) into the following 
system under the condition of steady flow (see appendix A) 


where 


dv 
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( 11 ) 

(12a) 

(12b) 

(12c) 

(12d) 

(12e) 

(12f) 
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Equations (9) to (11) contain five dependent variables ^p, p, a^, v^, and v^j. To 
solve for these variables, two additional equations are required. The integrated form of 
the energy equation for an ideal gas is used 

/ y \ V ^ + V , ^ 

^ ?- = H (13) 

P\y -1/2 

2 

along with an equation of state for an ideal gas relating a , p, and p 


a2 = VP 
P 


(14) 


Note that equations (13) and (14) are valid for ideal gases only and would have to be 
replaced by other appropriate equations if equilibrium chemistry were to be considered. 

Equation (9) is indeterminant along the stagnation streamline for the case of axi- 
symmetric flow because of the term (cot 0)v^; thus, it must be replaced by its limiting 
form along this line. (See appendix B.) For j = 1 and ^ = 0, 


where 


377 



3v , ^ 

2G, —2 + G 
, 3 a| 



(15) 


(16) 


The Method of Lines 

The method used to obtain solutions to the present blunt -body problem is patterned 
after the method presented in reference 19, for the solution of nonlinear conical flow 
problems, wherein it was called the method of lines (MOL). The application of the 
method of lines in the present problem proceeds as follows. The region of interest in the 
1,77 -plane (0 ^ ^ ^ 1, 0 ^ 77 ^ 1) is divided by NM straight lines parallel to the 77 -axis 
(see fig. 2), the line n = 1 being coincident with the stagnation streamline in the physical 
plane and the line n = NM being coincident with the line cp = in the physical 

plane. The boimdary 77 = 0 represents the body surface and the boundary 77 = 1 repre- 
sents the bow shock wave. This arrai^ement divides the region in NM - 1 strips with 
the NM lines forming the boundaries of the strips. In the present application, the 
spacing between lines is equal although this is not an inherent requirement of the method. 
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At each strip boundary or line, the system of partial differential equations (9) 
to (11) is reduced to a set of ordinary differential equations by replacing the ^-derivatives 
with finite -difference approximations. The resulting system of coupled ordinary differ- 
ential eqviations can then be integrated from an assumed bow shock wave to the body sur- 
face. This process must be repeated until a bow shock wave is found which yields the 
correct surface boundary conditions. The problem is formulated as a "direct free bound- 
ary problem" which has some features that are obviously similar to the inverse problem 
described in references 2 to 8. 

Finite -difference approximations.- Although several types of finite -difference 
approximations have been tried, one which approximately takes into account the physical 
zone of dependence in the flow has been found to yield much more satisfactory results 
than the other schemes that have been tried. In this scheme, at each step of the integra- 
tion, when the Mach number at a point on the nth line is subsonic (M < 1), a central differ- 
ence (ref. 21) is used to represent the ^-derivative. Thus, for M < 1, the derivative of 
some general quantity F becomes 


9F ^ ^n+1 ' ^n-1 
3 ^ ~ 2 


(17) 


When the Mach number is supersonic (M > 1), a three-point backward difference (ref. 21) 
is used; thus, for Mil 

3F ^ ^^n ' '^^n-l ^n-2 

9 ? 2 

Both expressions are formally second order accurate although the first expression has a 
smaller error term than the second. Using the more accurate central difference for- 
mula (eq. (16)) throughout the field has not been found as successful because it allows an 
improper direct coupling between downstream and upstream points in the supersonic 
region which slows convergence. 

Extension of MOL to three dimensions .- It is conceptually straightforward to extend 
the method of lines to treat three-dimensional blunt-body problems. For this case, lines 
would be constructed in a number of meridional planes around the body similar to those 
constructed in the ^,?]-plane described previously. This procedure would divide the flow 
field between the shock and the body into a number of regions or "cells" with the lines 
forming the edges of the cells. Along each line, the system of governing partial differ- 
ential equations would be reduced to a system of coupled ordinary differential equations by 
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replacing the partial derivatives in both the longitudinal and circumferential directions by 
finite -difference approximations similar to those described in the previous section. From 
this point, the similarity to the two-dimensional or axisymmetric problem should be 
obvious. The major advantage that the method of lines holds for three-dimensional prob- 
lems is that flow field data need not be stored at grid points between the shock and the 
body since a system of ordinary differential equations are integrated from the shock to 
the body to obtain the solution. This method leads to a reduction in computer storage 
requirements by approximately an order of magnitude compared with fully discrete 
methods (such as time -dependent methods) when the entire shock layer must be 
descretized. 

Integration of the ordinary differential equations . - The system of differential equa- 
tions (9) to (11) is integrated simultaneously along each line with the derivatives on the 
right-hand side being evaluated from equations (17) and (18). The finite -difference 
approximations cause the resulting differential equations at each line to be coupled to the 
differential equations on adjacent lines; thus, the system of differential equations becomes 
a set of 3NM coupled ordinary differential equations. Of the many methods available for 
integrating sets of ordinary differential equations (see ref. 22), a two-step predictor- 
corrector method proposed by Stetter has been used. Following this method, the pre- 
dictor step is 


pT7-AT]^„rj ^ 

n n 2 


and the corrector step is 




3(F’)J - 

- (1 - Vo)f (('F') 


rj-Ai) 

n 



(19a) 


(19b) 


where yq was taken as 0.65. This method is formally second-order accurate and very 
stable; however, it is a two-step method and requires some other starting solution to com- 
pute the first step. In the present study, a simple Euler-predictor trapezoidal corrector 
was used to start the integration and equal size steps were taken in the r; -direction 
throughout the entire integration cycle. Typically 5 to 10 integration steps were taken 
between the shock and the body. 
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Symmetry and Boundary Conditions 


The flow is symmetric about the stagnation streamline; thus, symmetry conditions 
are imposed along the line | = 0 (t] = 1) 


a,, 

V . = ^ = ^ 

0 34 04 


= 0 


( 20 ) 


At the start of each integration cycle, the shock jump relations are used to deter- 
mine all properties on the downstream side of the shock wave at ?7 = 1. (See appendix C.) 
The initial shock wave shape is assumed and then iterated until the surface boundary con- 
ditions are satisfied. The procedure used to correct the shock shape during the iteration 
process is described in a subsequent section. 

Along the body surface, the normal component of velocity must vanish. Thus, at 
the body surface (?] = 0) the following relation must hold; 


Vb = ^b - ^b = ° (21) 

The shock wave shape is initially assumed and iterated upon until this boundary condition 
is satisfied on each line. 

The downstream boundary (4 = 1) must be chosen so that at this boundary the flow 
is totally supersonic, and on this boundary three-point backward differencing is used. 

(See eq. (18).) 


, Determination of Shock Wave Shape 

The iteration procedure for determining the shock wave shape is the same as that 
used in reference 19 for conical flows. The shape of the assumed shock wave is given by 
the known numerical function r(7] = 1, 4 ) = r_(4); thus, this function can be differentiated 
numerically, using either equation (17) or (18), to determine drg^d4. Knowing r^(4) 
and drgy^d4 at each line allows the flow properties ^p, p, a2, Vj., and v^j on the 

downstream side of the bow shock wave to be determined from the shock jump relations. 
(See appendix C.) These flow properties can be used to evaluate the right-hand side of 
equations (9) to (11) to start the numerical integration procedure at 77 = 1 and proceed 
inward toward the body surface (?] = 0). Only the correct shock shape will satisfy the 
surface boundary conditions (eq. (21)). A corrective procedure must be used to drive 
the shock shape in a direction that will eventually satisfy these boundary conditions at the 
body surface. 
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Newton iteration for shock shape.- The shock wave is given by the equation 


rg(?) = + Ar(^) 


(22) 


Since the body shape r^(|) is always known (direct problem), the shock shape is 
expressed numerically by a set of NM shock-layer thicknesses Ar^^, one for each line. 
For each line, there is a normal component of velocity at the surface A Newton- 

type iteration procedure will be used for correcting the NM values of Arj^ so that 


max 




where e is a small prespecified convergence criterion (typically € = 10“^). The steps 
in the procedure are as follows: 

(1) Obtain an initial set of values for Arj^ either from experience or from a solu- 
tion to a previous problem. 

(2) Determine properties on the downstream side of the shock wave at each line. 

(3) Use the results of step (2) as initial conditions to start the numerical integration 
of the system of ordinary differential equations from the shock to the body surface and 
evaluate Vj^ at each line. 

(4) Determine maxj^Vjj)^|. If maxj^V| 3 )^| ^ e, the solution is considered con- 
verged; if not, proceed to step (5). 

(5) Successively perturb the value of the shock-standoff distance Ar^^ at each line 

to (1 + 6)Arjj, where 6 is a small number ^typically 6 = 10"® or 10"'^), and repeat 
steps (2) to (4). This procedure results in a small change in with each pertur- 

bation that can be related to a change in the shock-layer thickness at each line 


’(''b) 


n 




(n= 1,2,3,. . .,NM) 


After completion of the NM perturbations, a NM x NM matrix of "influence 
coefficients" is generated which can be used to correct the shock shape in a way so that 
(^b) driven toward zero. 

(6) Solve the first-order linear system 






(n = 1,2,3,. . .,NM) (23) 
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to obtain the corrections A^Arj^ to the value of the shock-layer thickness at each line 


which will tend to drive 


(''A 


toward zero. 


(7) By using the new shock shape 



new 



(24) 


begin a new cycle at step (2) and continue until the test at step (4) is satisfied. Note that 
each complete cycle requires one initial integration and NM perturbations (or NM + 1 
integrations). 

In reference 19, it was shown that a modified Newton iteration procedure could 
often be used to greatly speed up the convergence process. In this procedure, when 
max j(V^^ I ^ e (where e is a small number. Usually around 0.02 to 0.05), the pertur- 
bations were eliminated and the corrections to the shock-layer thickness at each line 
A^Ar^jj were evaluated from the matrix of old influence coefficients. This procedure 
reduces the number of integrations per cycle by NM and although it generally requires 
more cycles, it will usually require fewer total number of integrations (and thus less 
work) to obtain a converged solution. 

Initial shock shape .- The ability to obtain a converged solution for certain body 
shapes and flow conditions depends strongly upon the accuracy of the initial guess for the 
shock shape. However, for a sphere (or circular cylinder) and a paraboloid at moderately 
large Mach number > 5), converged solutions can usually be obtained by starting 
with a shock shape represented by a constant shock-layer thickness Ar^^ at each line. 

By using the solution thus obtained, it is usually easy to "build" the solution for other 
desired body shapes or flow conditions with the procedure outlined as follows. 

(1) A converged solution is obtained for a problem similar to the desired problem 
but differing in some parameter F where F could represent M^, y, or B^^. See 
section on body geometry. 

(2) A small change is made in F and a new solution is computed for 

Fnew = ^old initial value of AF is usually very small so that the old 

shock shape is a good approximation to the shock shape for the new case. 

(3) The changes in the shock-layer thickness Ar^j are obtained at each line 
between the new and old case and the derivative d Ar^j^dF is computed numerically. 

(4) A new set of approximate shock-layer thicknesses are obtained from the 
assumed linear relationship 
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+ 


(25) 



new 




AF 


and a solution is computed for a new value of the parameter F^^^^ = F^^^ + AF by using 
the shock shape obtained from equation (25) as an initial guess. 

(5) Steps (3) and (4) are then repeated until the desired problem has been solved. 
This procedure is easily adapted to automatic computation and is very efficient since it 
insures that the initial guess at the shock shape is a reasonably good approximation. 

This procedure has been built into the present computer program and is used extensively 
for obtaining solutions to the more difficult problems. 


Determination of Surface Properties 

From equations (9) to (11), it should be observed that if g or Gj becomes zero, 
a singularity occurs in the differential equations so that the numerators of equations (9) 
to (11) must also become zero if the 77 -derivatives are to remain finite. Numerically, 
evaluation of the tj - derivatives becomes very difficult when g or becomes small 
and the integration process must be stopped at this point. When the assumed shock shape 
is far from the converged shock shape, it may not be possible to integrate the differential 
equations all the way to the body surface without a singularity of the type mentioned above 
occurring. When this happens, the velocity components normal to the body 

V = Vr cos - v^ sin 

are extrapolated to the surface along each line by using a second-order accurate extrapo- 
lation technique which yields values of that can be used to correct the shock shape 

for successive iterations. As the shock shape approaches convergence, this type of singu- 
larity generally does not occur (if the last line does not extend too far into the supersonic 
region) and the integration of the differential equations can proceed to the surface. 

In evaluating flow properties on the surface, only the predictor step (eq. (19a)) in the 
numerical integration routine is used since the equations for the 17 -derivatives (eqs. (9) 
to (11)) are singular at | = 0. However, the predictor step is second order accurate and 
thus the results at the body surface have the same formal order of accuracy as the results 
away from the body surface. 


BODY GEOMETRY 

The body geometries that are considered in the present paper have a projection in 
the x,y plane that is described by the equation for a general conic section 
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(26) 


y2 = 2Rj^x - B^jx2 


where is the radius of curvature at x = 0 and is a bluntness parameter. This 
equation can be placed in a more convenient form if x and y are nondimensionalized 
by dividing by the radius of curvature R^j. Thus, 


y2 = 2x - B^x2 


(27) 


The bluntness parameter B^ characterizes the eccentricity of the conic section. Its 
significance is better understood if it is noted that < 0 generates a hyperbola, 

Bjj = 0 generates a parabola, and B^, > 0 generates an ellipse, with B^ = 1 for the 
special case of a circle. For an ellipse, B^ is related to the ratio of major to minor 
axes (b/a) by the relation 



An axisymmetric body is generated by revolving the curve described by equation (27) 
about the x-axis. By using the following equations to define x and y (see fig. 3), 



(29) 


an equation for a general conic in polar coordinates is obtained 


+ ~ 4AC 

2A 


(30) 


where 


2 9 

A = sin 0 + B^ cos‘^0 
B = 2 cos cf)(l - Bfe^p) 
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For the results presented in this paper, when > 0, Xp is taken as 

(31) 


which places the pole of the coordinate system at the point of intersection of the major 
and minor axes of an ellipse. (See fig. 3.) When ^ 0, Xp is set equal to 1 which 
places the pole of the coordinate system at the center of curvature of the nose. 

Thus, equation (30) is very convenient for describing the body shape since a wide 
variety of bodies can be obtained by varying a single parameter B^. 


RESULTS AND DISCUSSION 

In this section, results of computations by the method of lines (MOL) for several 
smooth body shapes are compared with experimental data and the results of other com- 
putational techniques to demonstrate the accuracy and utility of the present method. 


Sphere 

Pressure distributions over a sphere for y = 1.4 are presented in figure 4 for 
free-stream Mach numbers of 2, 3, 4, 6.05, and 8.06. For each case 9 lines were 
used and 10 equally sized integration steps were taken between the shock and the body. 
Experimental data from references 23 and 24 and computed results from reference 5 are 
presented for comparison. The pressures computed by the MOL agree well with the 
experimental data for each Mach number considered. The computed results from refer- 
ence 5 (using an inverse method) agree with the other data in the stagnation region, but 
underpredict the pressure as the sonic line on the body is approached (0 = 41° to 48°), 
particularly at low Mach numbers. 

The shock-layer thicknesses for these cases are compared in figure 5, where it is 
seen that for M^ > 3, the shock -layer thickness predicted by the MOL agrees closely 
with those from the other data. For M^ = 2, the MOL predicts a slightly thicker shock 
layer (by as much as 6 percent when compared with the experimental data of ref. 23); 
however, the predicted surface pressures were found to agree with the experimental data. 
(See fig. 4(a).) Since the shock-layer thickness changes rapidly at low Mach numbers, 
this behavior could easily be due to a small change in the effective Mach number for the 
tests which would change the shock-layer thickness but have little effect on the measured 
pressure distribution. 

Experience gained through application of the MOL to flow over a sphere and other 
body shapes indicates that for y = 1.4, M^ = 2 represents a practical lower limit for 
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using the method. Although solutions can be obtained at slightly lower Mach numbers, 
the computations become tedious and require more "art than science" to be successful. 

In references 19 and 20 it was pointed out that when the MOL is applied to Laplace's 
equation, roundoff errors grow exponentially like 10"! exp(NX) where j is the number 
of decimal figures in a computation, N is the number of lines, and X is distance along 
a line from the initial data line. Since the physical shock-layer thickness grows rapidly 
with decreasing Mach number (for Mach numbers approaching 1), it should be expected 
that this inaccuracy would lead to convergence problems as the Mach number is reduced. 
This has been found to be the case in the present study. 

For the = 8.06 condition, a matrix of solutions have been computed by using 
5, 7, 9, or 11 lines; 5 or 10 integration steps between the shock and the body; and a full 
or modified Newton iteration scheme. The results obtained from these computations are 
compared in tables I, II, and HI. Table I lists the conditions for the 16 cases and some 
results including number of global iterations required for convergence, total number of 
integrations, and computational time on a CDC CYBER 175 computer for each case. 

Table n lists the pressures on the body for the 0 locations associated with the solution 
for nine lines and table III lists the shock-layer thickness for these same 0 locations. 

In tables II and III, results are presented only for the full Newton iteration scheme since 
there is almost no difference between those and the results for the modified iteration 
scheme. All solutions were for ~ were started with a constant shock- 

layer thickness Ar = 0.140 at each line (concentric shock). 

Note first that the results obtained with the modified iteration scheme (e > e, 
cases 9 to lo) in general required more global iterations but fewer total integrations than 
with the full Newton iteration scheme (e = e, cases 1 to s). Since the actual amount of 
work involved in obtaining a solution is directly proportional to the total number of inte- 
grations, the modified iteration scheme required less computer time. However, the 
reduction in computer time is relatively small (of the order of 10 percent or less) and, 
in general, has not been found to be large enough to offset the disadvantage that for some 
problems the modified iteration scheme does not converge. Thus, the full Newton itera- 
tion scheme has been used exclusively in the remainder of this paper. 

With 11 lines, when 10 steps were taken between the shock and the body (case 1), 
the solution did not converge. When 5 steps were taken (case 5), the solution passed the 
convergence criterion, but the entropy on the surface was not constant near the last line 
and indicated a poor solution. The pressure on the surface at the last line is too high 
(table II) even though the shock-layer thickness is approximately correct (table HI). As 
discussed previously, roundoff errors grow exponentially with the number of lines so that 
as the number of lines is increased, it becomes more difficult to obtain a good converged 
solution. From experience, it has been found that 9 is a practical upper limit for the 
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number of lines. Some solutions have been obtained with more lines, but the computations 
become tedious. For all cases, the last line cannot extend too far into the supersonic 
region, but when the number of lines becomes large, the placement of the last line rela- 
tive to the sonic line (which is not known a priori) becomes critical. 

By comparing the results in tables n and m, it is seen that the solutions to the 
present problem with as few as 5 lines and 5 steps between the shock and the body are 
very accurate. In general, it has been found that using up to 9 lines with 5 to 10 integra- 
tion steps between the shock and the body provides very accurate results for most smooth 
blunt -body problems considered here. 

In figure 6, constant-density lines are presented for flow over a sphere in nitrogen 
with = 5.017 along with experimental data from reference 25 and computed results 
from reference 7. The results computed by MOL are in excellent agreement with those 
of reference 7 and reasonably good agreement with the experimental data except in the 
stagnation region where because of the uniformity of the flow, accurate measurements 
are very difficult to obtain. 

In figure 7, constant -density lines are presented for flow over a sphere in argon at 
= 5.329 along with experimental data from reference 25 and computed results from 
reference 7. There is excellent agreement with the computed results of reference 7, but 
the agreement with the experimental data is poor. The reason for the disagreement with 
the experimental data is not known; however, the good agreement with the computed 
results of reference 7 supports the accuracy of the present method. 

The MOL was one of several computational methods used in reference 26 to com- 
pare with experimental shock shapes over a sphere covering a density ratio across the 
normal shock from 4 to 19. These density ratios represent a range of effective y's 
from 1.667 to 1.067. In that report, it was shown that the results obtained by the MOL 
agreed well with the experimental shock shapes and with most of the shock shapes pre- 
dicted by the other computational methods. 

Paraboloid 

Solutions for the pressure distribution over a paraboloid with y = 1.4 are presented 
in figm-e 8 for Mach numbers of 10 and 3. Computed results from reference 5 are included 
for comparison. The solutions agree very closely except in the vicinity of the outflow 
boundary where the two solutions deviate slightly. These are typical for the com- 

parisons presented previously for a sphere (see fig. 4) except that for a sphere, experi- 
mental data were also available which substantiated the MOL results. 

The variations of shock-layer thickness with <p for these cases are compared in 
figure 9. The shock -layer thicknesses for the two solutions agree very closely. 
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Ellipsoid 


Pressure distributions for an ellipsoid with b/a = 0.5 and y = 1.4 are presented 
in figvtre 10 for Mach numbers of 8.06, 6.05, and 4. Experimental data from reference 23 
are included for comparison. The computed results agree very closely with the experi- 
mental data for each Mach number. 

Pressure distributions for an ellipsoid with b/a = 1.5 and y = 1.4 are presented 
in figure 11 for Mach numbers of 8.06, 6.05, 4.0, and 3.0 along with experimental data 
from reference 23. The computed pressure distributions agree closely with the experi- 
mental data. The variations of shock-layer thickness with 0 for this body shape are 
presented in figure 12 for Mach numbers of 8.06 and 3.0 where it is seen that the MOL 
results agree closely with the experimental data. 

In figure 13, the shock standoff distance at the stagnation point is presented for a 
body bluntness parameter = (b/a)^ ranging from 0 (paraboloid) to 9 (3:1 ellipsoid) 
for M_^ = 10, 6, and 3. The only data available for comparison are the computed results 
from reference 7 which cover a limited range of Bjj and M^. The shock standoff 
distances computed by the MOL compare well with the available data from reference 7. 
This figure illustrates the wide range of body shapes for which the MOL can be applied. 

Circular Cylinder 

Results for the two-dimensional flow over a circular cylinder are presented in 
figures 14 and 15 for Mach numbers of 8.06 and 3.0. Figure 14 presents MOL surface 
pressure distributions compared with experimental data from reference 23 and com- 
puted results from reference 6 (inverse method) and from a perfect gas version of the 
method presented in reference 27 (time-dependent method). The MOL results are in 
good agreement with those of the other data. 

Figure 15 presents MOL shock-layer thicknesses compared with the computed 
results from reference 6 and the method of reference 27. Experimental data for the 
shock-layer thicknesses for this case are not available. All results agree well at 
M^ = 8.06. At M^ = 3, there are much larger differences among the results. The 
results computed by MOL and the method of reference 27 are in reasonably good agree- 
ment, but the shock- layer thicknesses from reference 6 are much smaller. Whatever 
the reason for the differences in the shock-layer thicknesses, the pressures on the sur- 
face are unaffected. (See fig. 14(b).) 

CONCLUDING REMARKS 

The method of lines (MOL) has been used to obtain solutions to axisymmetric and 
two-dimensional inviscid flow over smooth blunt bodies. Comparisons with experimental 
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data and the results of other computational methods have demonstrated that very accurate 
solutions can be obtained by using relatively few lines with this approach. The method 
is efficient; typical converged solutions have been obtained in 100 to 150 total integrations. 
The method of lines is semidiscrete and has relatively low core storage requirements as 
compared with fully discrete methods (such as time -dependent techniques) since very 
little data are stored across the shock layer. This latter feature is very attractive for 
three-dimensional problems because it enables core storage requirements to be reduced 
by approximately an order of magnitude. 

The disadvantage of the method of lines is that roundoff errors increase exponen- 
tially with number of lines; thus, the total number of lines that can be used for a particu- 
lar problem is limited. In the present study it was found that 9 lines was a practical 
upper limit for two-dimensional and axisymmetric problems. This limits application of 
the method to smooth body geometries where relatively few lines would be adequate to 
describe changes in the flow variables around the body. 

Extension of the method to three dimensions is conceptually straightforward; how- 
ever, three-dimensional application would also be limited to smooth body geometries 
although not necessarily to a total of 9 lines. 

Langley Research Center 
National Aeronautics and Space Administration 
Hampton, VA 23665 
January 27, 1978 
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APPENDIX A 


TRANSFORMATION OF FLOW FIELD EQUATIONS 

Consider the transformation of the flow field equations from physical space (r,0) 
to the computational space ( 77 , |). The flow field equations in physical space (eqs. (2), 
(3), and (5)) can be written in the following expanded form for steady flow: 


r 9r r 90 r p dr 


(Al) 


r 8r r 80 r pr 80 


(A2) 


r 8r r 80 '^\8r r 80 r r ■’ r 


(A3) 


Using the transformation operators defined by equation (7), equation (Al) transforms as 
follows: 






r\Ar 977 / r V '0 877 8 | / 


'A 


ifiaUo 

P\Ar 977/ 


(A4) 


which after solving for Qv^jdri yields equation (10) 


!!l= l( 1 9p , ^ 0^0 ^^r ^0 

877 g\Arp dr] ^ r 8 ^ r 


(A5) 


Similarly, equation (A2) transforms as follows: 


r 0 


pr 




(A6) 
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which after solving for 



yields equation (11) 


9?7 "g 


If. 

pr\^<l) <t> &r]} \ T j r 


(A7) 


Finally, equation (A3) transforms as follows: 


V JJL M 2E a2p 

^\At drij r \'<}> drq ^0 


1 9v^ 1 / av , 9v A 

— + ^ ^ 

Ar 9rj r\'0 977 ^0 9| ^ 


1+J.V jCOli 

r r •’ r 0 


= 0 


(A8) 


Now substituting these expressions for and bv ^j&r] and solving for 9p/9rj 

yields equation (9) 


iP= - ifog^ + 03^ + 0.— + Gc 
977 Gj \ 2 3| 3 4 9| 5J 


Note that g, Gj^, G 2 , Gg, G^, and Gg are given by equations (12). 
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APPENDIX B 


LIMITING FORM OF FLOW FIELD EQUATIONS 
ON STAGNATION LINE 

For axisymmetric flow equation (9) 


dr] 



3v 


-^ + G. 



is indeterminant along the stagnation streamline. Thus, a limiting form of this equation 
must be derived. Symmetry conditions apply along the stagnation streamline; thus, 


4 > 0 


(^ = 0 ) 


Applying these conditions to this equation yields 




where 




(Bl) 


(B2) 


(B3) 


and 


G 


5 ■ 



+ g cot 



(B4) 


The term cot (f> in the expression for Gg is indeterminant along the stagnation 
streamline = 0); thus, the limiting form of this expression must be used 
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av , 3v , 
lim cot 0 V , = — £ = I , — ^ 
0^0 ^ (P d(p ^<t> 


(B5) 


By using this result, equation (B4) becomes 




= G, -^ + Gc 
3 3| E 


(B6) 


where 


G5 = i-£ 

D rg 


(2v,g) 


(B7) 


Substituting equation (B6) into equation (Bl) yields equation (15) 


av 


2G, + G 


dv 


9? 


(4 = 0) 


which is determinate along the stagnation streamline for axisymmetric flow. 
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APPENDIX C 


SHOCK JUMP CONDITIONS 

Consider the flow at the bow shock wave depicted in sketch (a). 



The shock wave angle B is given by the equation: 

s 




(Cl) 


where 




The polar velocity components at the shock wave are given by the equations 

W =‘K) cos cos sin Cg 

s s 

(^4>) = (^^n) si" ?s + ^00 cos cos Cg 

s s 

The oblique shock wave relations can be written in the form (ref. 28) 
p^V sin B = p (v ) 

f^oo oo Hg jj) 


(C2) 


(C3) 

(C4) 


(C5) 
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+ P„V J + p^(v J 


(’n): 


I h , ’ 

h + = + 

oo 2 s 2 


The nondimensional free-stream variables have the following values: 


V ^ P =1 

* oo ^°o 




^oo « - 2 

y^«r 


(C6) 


(C7) 


(C8) 


and for a perfect gas the density on the downstream side of the shock wave is given by 
the following equation: 


(y H- 1)M^^ sin^^g 
' (y - 1)M^2 + 2 


(C9) 


Thus, knowing y and and the shock wave angle which can be determined 

from equations (Cl) and (C2), the properties on the downstream side of the shock 
wave ^Vjj^ , Pg, and h^ can be calculated from equations (C5), (C6), and {Cl). 

Finally, ^v^.^ and ^v^^ are calculated from equations (C3) and (C4), 
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TABLE I.- CONDITIONS FOR TEST CASES 


[ 

Moc = 

8.06 

; y 

= 1.4 



Case 

NM 

NET A 

€ 

€ 

Iterations 

Integrations 

Time, 

sec 

1 

11 

10 


0.001 

1 

No con 

vergence 

— 

2 

9 







17 

170 

7.123 

3 

7 







16 

128 

4.349 

4 

5 

10 





13 

78 

2.031 

*5 

11 


5 





8 

96 

2.623 

6 

9 







10 

100 

2.268 



, 






10 

80 

1.486 




5 



.001 

10 

60 

.862 


11 

1 

0 



.05 

No convergence 

— 

10 

9 







21 

156 

6.564 

11 








21 

126 

4.313 

12 


10 





17 

72 

1.866 

*13 

11 


5 





11 

77 

2.239 

14 

9 







13 

76 

2.011 

15 








15 

71 

1.402 

16 








13 

53 

.805 


*Surface entropy check failed on last line. 
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TABLE n.- SURFACE PRESSURE DISTRIBUTIONS 


FOR TEST CASES 


0, 

deg 

Surface pressure distributions for case - 

1 

2 

3 

4 

5 

6 

7 

8 

0 

5.625 

11.250 

16.875 

22.500 

28.125 

33.750 

39.375 

45.000 

— 

0.9246 

.9141 

.8831 

.8330 

.7676 

.6894 

.6037 

.5148 

.4278 

0.9246 

0.9246 

0.9241 

0.9241 

.9136 

.8826 

.8324 

.7673 

.6891 

.6025 

.5132 

.4261 

0.9242 

0.9243 


.8832 



.8829 




.7674 

.7690 

.7669 

.7671 

.7687 


.6013 



.6012 




.4274 

.4269 

.4397 

.4271 

.4272 


TABLE m.- SHOCK- LAYER THICKNESSES FOR TEST CASES 


deg 

Shock- layer thicknesses for case - 

1 

2 

3 

4 

5 

6 

7 

8 

0 

5.625 

11.250 

16.875 

22.500 

28.125 

33.750 

39.375 

45.000 

— 

0.1401 

.1409 

.1433 

.1473 

.1532 

.1614 

.1723 

.1863 

.2040 

0.1402 

0.1406 

0.1407 

0.1407 

.1415 

.1440 

.1480 

.1539 

.1622 

.1731 

.1871 

.2048 

0.1409 

0.1414 


.1439 



.1447 




.1533 

.1537 

.1540 

.1541 

.1545 


.1730 



.1740 




.2043 

.2055 

.2048 

.2053 

.2067 
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